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1. Introduction 

As the latest generation of gravitational wave detectors becomes operational, the 
problem of faithfully simulating the evolution of binary systems of compact objects, 
black holes in particular, has become increasingly important. While post-Newtonian 
(PN) approximations can be used in the first stages of the life of a binary black hole 
(BBH), when the two objects get close and are rapidly orbiting around each other 
only solutions to the full non-linear Einstein equations can provide the desired level of 
precision. Due to the complexity of such equations, these solutions can only be achieved 
by means of numerical algorithms. These results are of particular interest for laser- 
interferometric observatories since BBH will be highly relativistic when entering the 
sensitivity range of the detectors. 

Such simulations pose a hard and challenging problem. Until recently they tended 
to fail after a very short time due to instabilities (25] which resulted in exponentially 
growing run-away solutions. Fortunately, tremendous progress has been achieved over 
the last two years [HI [32l H31 [3l HH [H [2H [151 EH HSl El [33l H2l [221 S21- Within the 
moving puncture approach and also with the generalized harmonic system, it is now 
possible to evolve BBH systems through several orbits and the subsequent merger and 
ringdown phases. 

Stable and accurate, modern BBH simulations require large computer resources 
and even modest size runs are performed on supercomputers involving dozens or even 
hundreds of processors. The goal of this paper is to showcase the ability of the 
code BAM [12] to evolve certain BBH systems on workstations, providing results 
of comparable quality to those obtained in simulations using much larger computer 
systems. Workstations with similar characteristics to the ones used here are reasonably 
affordable (less than $3000 USD at the time of publication). BAM provides grid 
structures composed of boxes of increasing resolution near the center of the grid. In 
the case of binaries, the highest resolution boxes are placed around each black hole. 
The boxes track the holes in their orbits until the final merger, when a single set of 
levels surround the black hole remnant. A direct consequence of this grid structure is 
the efficient use of computational resources as it will be detailed in section [31 BAM 
currently handles fourth order accurate evolutions. 

The end result of a BBH merger is a larger black hole. This final black hole could 
in principle be non-spinning, however the conditions for this to occur are very unlikely: 
any astrophysically realistic scenario would lead to a spinning object. We test BAM by 
simulating a BBH with identical black holes with intrinsic spins S/m = 0.75 parallel 
to the orbital angular momentum. We choose high-spin binaries since they require 
very high resolution near the black holes which is currently difficult to achieve for most 
numerical codes of this type. The time evolution of the BBH system is achieved through 
the moving punctures method [131 Ell- 
in order to test the quality and accuracy of the simulations, we concentrate on 
the measurement of the final black hole mass and angular momentum. To do that, 
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Table 1. Initial data parameters. Here nib is the bare mass parameter of each puncture 
and M — 2m is the sum of the ADM masses m measured at each puncture. The holes 
have coordinate separation D, with puncture locations (0,±-D/2,0), linear momenta 
(=FP,0,0), and spins (0,0,5) with S/m 2 = 0.75. We also list the 2PN estimates for 
the ADM mass M^ DM , the ADM angular momentum J^ DM and the angular velocity 
O. These quantities are shown using two different scaling factors (M and M^ DM ) for 
easier comparison with work done by other groups. 



we implemented algorithms based on the conversion of surface integrals (at the core 
of the definition of the ADM mass and angular momentum) to volume integrals using 
Gauss' theorem. These calculations are studied and compared with alternative ways of 
measuring these global quantities. 

Section [2] presents a brief description of the equations and the initial data sets and 
describes the details of BAM's numerical grid structure [12] and the algorithms used 
to calculate the mass and angular momentum. Sections [3] and |4] presents performance 
and convergence tests respectively. Section 14.21 compares alternative calculations of the 
mass and angular momentum and our results are discussed in section [SI 



2. Evolution using the moving punctures method 

2.1. Initial data 

In order to start our simulations we need initial data for spinning BBH with equal masses 
and spins. Since we will employ the moving punctures approach in our evolutions we 
will use standard puncture initial data [10] with the momentum and spin parameters 
in the extrinsic curvature given by 2PN estimates [26J. It is sufficient to use 2PN 
estimates because standard puncture data are inconsistent with PN theory beyond 
(v/c) 3 [3U HH US]. These parameters along with 2PN estimates for ADM mass M ADM , 
ADM angular momentum J ADM and angular velocity Q are shown in Table [U 

The coordinate distance D and the momentum and spin parameters P and S 
directly enter the Bowen-York extrinsic curvature, while the bare mass parameter 
is obtained from the condition that the ADM masses measured at each puncture 
should be m = M/2. This implies that each black hole has an individual spin of 
m? = (m%) 2 = 0-75, where as in [39], H0[ [2] we assume that m is a good approximation 
for the initial individual black hole masses. Note that these data are very close to the 
values used by Campanelli et al. [IS] • If we express everything in terms of the PN ADM 
mass, the largest difference is that our bare mass parameter is about 1% lower. 

To complete the definition of the initial data, we also need to specify initial values 
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for the lapse a and shift vector (3 l . At time t = we use 

V ri r 2 / 
/3 l = 0, 

where ta is the coordinate distance from puncture A. Both lapse and shift are updated 
by evolution equations depending on the physical variables, as described below. 

2.2. Evolution of gravitational fields 

We evolve the initial data with the BSSN system [35J E]. In the case of BSSN, the 
3-metric ^ is written as 

9ij ^ ^9ij 

where the conformal metric has unit determinant. In addition, the extra variable 
T = -dfg ij 

is introduced where g^ is the inverse of the conformal metric. Furthermore, the extrinsic 
curvature is split into its trace free part and its trace K, and given by 

These variables are evolved using 

d (p = - -aK, 
o 

dogij = - 2aA ij , 
d A l3 = e-^l-DiDja + aR i3 ] TF 
+ a(KA ij -2A ik A* j ), 

d K = - D l D t a + aiAjAv + ^K 2 ), 



2 -f 
3 

2 



where do = d t — £p, Di is the covariant derivative with respect to the conformal metric 
c/ij, and "TF" denotes the trace-free part of the expression with respect to the physical 
metric, Xf- F = Xij — ^gijX^. The Ricci tensor R^j is given by 



Rij -^ij 

Rij = ~ \g lm did m ~ gi] + g k (idj)t k + f k f (ij)k + 



g lm \ 2T f(jf j)j TO + T k m T Hj j . 
- 2DiD j( f) - 
Ag ij D k <t)D k 4). 



R% = - 2DiD j( j) - 2g ij D k D k (j) + 4D, 
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The Lie derivatives of the tensor densities 0, and Ay (with weights 1/6, —2/3 and 
-2/3) are 



C P (j> =p k d k <p+-d k (3\ 

As&i = 9ijd k gij + g ik dj(3 k + g jk di(3 k - -gijd k (3 k , 
_ _, 2 ~ 

= Ay^y + AifcS,-/?* + Aj-fcdf/3* - -Aijd k fi k . 



As in [%Tl [12] we evolve the BSSN system as a partially constrained scheme, where 
both the algebraic constraints det(gr) = 1 and Tr (Ay) = are enforced at every 
intermediate time step of the evolution scheme. In addition, we also impose the first- 
order differential constraint P = —djCj 1 ^ by replacing all undifferentiated occurrences of 
P by —djg % i instead of using the evolved variable P. 

Note that for puncture initial data the BSSN variable has a divergence of the form 
logr^ at each puncture. Since a logarithmic divergence is relatively weak, the moving 
puncture approach consists of simply ignoring this divergence by putting it between grid 
points at the initial time. One option is to simply evolve the resulting initial data using 
a finite differencing scheme, which effectively smooths out any divergences, obviating 
the need for any special treatment of the punctures. This is the approach we have 
followed in our so called P-runs. Another option is to replace the BSSN variable by 
a new variable [13J 

X = e-'t 

which initially goes like r\ at puncture A. We use this second option in our C-runs. 

The second ingredient in the moving-puncture method is a modification to the 
gauge choice. We use a "1+log" lapse of the form [13] 



with 7] = 1.0/M for the P-runs. For the C-runs we used the modified gamma-freezing 
condition 



where advection terms have been added to all time derivatives, and where we choose 
r] = 2.0/M. 

2.3. Measurement of the mass and angular momentum 

The black hole resulting from a BBH merger is defined by its mass and angular 
momentum, thus an accurate measurement of such global gauge-independent quantities 
becomes critical. One way to estimate these quantities is evaluating the ADM mass and 



(d t - p { di)a = -2aK. 
For the shift, we use the gamma- freezing condition [Tl [T3] 




(1) 



(dt - P k d k )f3* = *-B\ (dt - (3 k d k )B l = (d t - P k d k )r - V B* 



(2) 
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angular momentum after the merger. They are defined as surface integrals on a surface 
arbitrarily far from the system [301 E] 

M ADM = 1 / ^ _ q j g nm g lr ^ (3) 
l07T Joo 

jf DM =^e i rf oo X l AldS n , (4) 

where e^- is the Levi-Civita tensor and (iS*™ = \sfg £nim dx l dx m . The estimation of the 
final black hole parameters in numerical codes is currently done in several different ways: 
1) by evaluating Eqs. (J3], [4]) as far as the grid size permits, 2) by measuring properties of 
the apparent, event or isolated horizons (see for instance [SUE]), 3) by estimating the 
amount of emitted energy and angular momentum in the form of gravitational radiation 
or 4) by converting the surface integrals (El Hj) to volume integrals using Gauss' theorem. 
This last method, which has been successfully employed in accretion disks around black 
holes (see for instance [20], [21]), binary neutron star systems [36] [18] [27] [28] [29] and 
single black hole spacetimes [43], is employed here and compared with results obtained 
from some of the other techniques. Some of the advantages of this method are 1) the 
reduction of the influence of noise generated at the outer boundaries, 2) the reduction of 
gauge drift effects, 3) it provides a real-time quality control factor at all times during the 
simulation and 4) it complements, and sometimes improves, the accuracy of alternative 
measurements. These advantages will be clarified in the following section with examples 
and comparison with some of the alternative methods. All the formulas used here are 
in Cartesian coordinates. The derivations in this section follow those of Yo et al. [43] 
and Duez [T9] . 

Since we are interested in applying these equations in formalisms that perform a 
conformal decomposition of the physical metric = e 4 * it is useful to transform 
them accordingly. Following [43], Eqn. ([3]) can be re- written as 

M ADM = 1 / [e * (Q- mr _ Q m - ir) + 4 ^ - m _ Q me <t>g lr ^ 
l07T Joo 

g nm 9 lr dS n 

= 77- ( f [9 lr (dig m r - d m9lr ) - 8 d m e*} g nm dS n 

1D7T Joo 

= ^ £(r n - Y l \ - 8 7D V) dS n , (5) 

with dS n = \\fg e n im dx l dx m , Y l - k the conformal affine connections, T l = —djg 1 ^ and D n 
the covariant derivative with respect to the conformal spatial metric. Similarly, Eqn. 
(T4]) becomes 

df DM = ^ ej" / x l A« m dS n . (6) 

07T Joo 

Note that we have not assumed that ^/g = 1, as is the case in the BSSN formalism. 

The method studied in this paper converts straightforwardly these surface integrals 
using Gauss' law with the only provision of excluding the parts of the grid immediately 
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surrounding the black holes. For an arbitrary vector field /j, Gauss' law adopts the form 
If n dS n =[ d n (y/g f n )dx 3 + Y,{ f n dS n , (7) 

where the first term represents the volume integral over all space except the parts 
enclosed by the closed surfaces dQk that surround each one of the k black holes. 
Numerical simulations like ours are performed using grids that cover a finite spatial 
volume V. Because of that, we can only provide estimates to the quantities defined in 
Eqs. (0 ED- We will call the calculations of the mass and angular momentum performed 
in our finite grid volumes My and Jy respectively. After the merger, the gravitational 
fields settle down in the Kerr geometry corresponding to the final black hole. Once 
this occurs, the values of My and Jy should approximate the corresponding mass and 
angular momentum of the black hole remnant. The convergence of this approximation 
is studied in the next section. 

In the case of the mass formula (jSJ), the direct application of (jTj) plus some algebra 
leads to 

M v = 777" / ( R + f " Tl nl ~ T lnm T nlm - 8 /3V) ^ d 3 x 

l07T JV 

+ 717-E/ (T n -T n \-8D n e*)dS n , (8) 



167T ^ Jdn k 

where R is the spatial Ricci scalar, D 2 = g mn D m D n , and the infinite volume Voo has 
been replaced by the finite volume V. The final step in the derivation of the mass 
formula is the use of the Hamiltonian constraint 

D^e* = -R+ —K 2 - — A mn A mn , 
8 12 8 

to eliminate in Eqn. (jSJ) the term proportional to D 2 e^ 

M v = -L f [(1 - e*) R + T n Y l nl - Y lnm Y nlm + 

l07T JV 

e 5 *( A mn A mn - ^K 2 )} ^ d 3 x 

+ 77^-E/ (r-fVsnV)^. (9) 

lbiT ~ Jan k 

In a similar manner, Eqn. ([6]) is converted to 



\- el I V H A l m + x l D n (e^A n J - x l A ns d m g r 
m Jv l 2 



d 3 x + 



JiV — 7— e.- 

87T 

7- e « m E£ x l AldS n . (10) 
8tt V Jdci k 



The momentum constraint 

D n (e^A n J = 2 - e^D m K 



is now used in Eqn. ( TTUI) . resulting in 



JlV ~M €u Jv" 



~i 2 - 1 - 

A m g x D m K — —x A ns d m g 



8tt €il ^£n k 



d 6 x + 

6 'I J 

x l A n m dS n . (11) 
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Figure 1. Schematic diagram of BAM's grid structure. 

Our calculations were based on the BSSN formalism for which the conformal 
transformation = g~ij is such that <\> = ln(^ 1 / 12 ). This condition imposes the 
algebraic constraints g = 1 and T l nl = which can be used to simplify even more Eqns. 
(Q and (fTTl) . These equations are the formulas used in this article with the "bar" fields 
replaced by the corresponding BSSN ("tilde") counterparts. 

Figure Q] shows the schematics of a typical grid structure. The darker shading 
indicates higher resolution. In order to avoid the coordinate singularity in the calculation 
of My and Jy, we exclude the region surrounding the black holes. Given BAM's mesh 
structure, it is particularly simple to choose the outer edge of one of the moving boxes 
(the innermost in Fig. [1]) as the boundary dQk- In the next section we provide results 
for different choices of these boundaries. Both numerical integrations (volume and 
surface) are performed using an extended version of the trapezoidal rule with fourth 
order convergence [31] . 

3. Code performance 

The numerical results discussed in this paper were obtained with the BAM code [TTl[T2] . 
This code is based on a method of lines approach using fourth order finite differencing 
in space and explicit fourth order Runge-Kutta (RK) time stepping. For efficiency, 
Berger-Oliger type mesh refinement is used [8]. The numerical domain is represented 
by a hierarchy of nested Cartesian boxes. The hierarchy consists of L + 1 levels of 
refinement, indexed by I — 0, . . . , L. A refinement level consists of one or two Cartesian 
boxes with a constant grid-spacing hi = ho/2 on level I. We have used here L = 9 
to 11 for the number of refinement levels, with the levels through 5 each consisting 
of a single fixed box centered on the origin (the center of mass). On each of the finer 
levels 6 through L, we initially use two sets of moving boxes centered on each black 
hole. When the black holes get close enough that two of these boxes start touching, 
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they are replaced by a single box. The position of each hole is tracked by integrating 
the shift vector. We have used this same set up with different resolutions to perform 
convergence tests. The notation used to describe these grid setups is as follows: the CI 
run is represented by 

• CI: Xr?=2[5 x 40 : 6 x 80][h w = M/56.9 : OB = 729M] 

where x represents the use of that dynamical variable, rj is the free parameter in 
the shift vector formula, 5 x 40 indicates that we have 5 levels with moving boxes 
of 40 x 40/2 x 40/2 points and 6 levels with non-moving boxes of 80 x 80/2 x 80/2 
points (the divisions by 2 are due to using quadrant symmetry). The resolution at the 
finest level (I = 10 in this case) is h w = M/56.9 and the outer boundary is placed at 
~ 729M from the origin. 

Finally, we note that BAM is MPI parallelized. When iV processors are used, 
each box on each refinement level is divided into N equally sized sub-boxes with added 
ghostzones. Each of these sub-boxes is then owned and evolved by one processor. The 
ghostzones are synchronized in the usual way after each evolution step. In this way, 
each processor owns exactly one sub-box of every mesh refinement box, which optimizes 
load balancing since then each processor works on the same number of grid points. For 
additional details about the version of the BAM code used here, see [T2] . 

We tested the BAM code by running simulations of a high-spin black hole binary 
system with individual spins S/m 2 = 0.75 aligned with the orbital angular momentum. 
We choose x—y as the binary's orbital plane which leaves the z component of the angular 
momentum as the only non-zero component of Jy. Table [2] shows the characteristics 
of the simulations performed for this article. The runs are grouped in those using \ 
(C-runs) as the dynamical variable and those using (P-runs). For the former (latter), 
we chose a value for the shift parameter rj of 2.0/ M (1.0/M). In general, the simulations 
performed here have higher maximum spatial resolutions than those in [12] , showing the 
moving punctures method's stability even when the coordinate separation between the 
grid points and the punctures is as small as M/320 (run P6) fl Another difference is 
that the grids used for the C-runs are larger than those of [12]. These extensions of the 
grid size and high resolutions did not excite any undesirable instabilities- in the case of 
C2, the simulation was run for the equivalent of over 20 orbital periods 

One of the strongest characteristics of the BAM code is its ability to perform good 
quality simulations with modest computer resources. Table [3] shows typical running 
times and memory requirements for the simulations of Table [2j Note that BAM performs 
faster after the merger, when only one set of non-moving boxes remains. All the 
simulations presented in this paper were run on dual processor workstations, namely 
a AMD Dual Opteron 2.2GHz workstation with 8Gb of memory and a Intel Dual Xeon 
2.6GHz workstation with 16Gb. The former computer can be purchased at the time of 
publication for less than $3000 USD. Note that, while none of the runs presented here 

| See [23j [7] for details on the spacetime geometry close to the punctures. 

§ Assuming a nominal orbital length of 114M, obtained from the initial data set angular velocity. 
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Table 2. Grid structure of the x (C) ancl (P) an( l runs. The x (0) runs used 
r/ = 2.0/M (77 = 1.0/M). 



required significantly more than 8Gb of memory, the cost of expanding the workstation 
memory capabilities has dropped considerably in the past year. Currently, a 1Gb 
memory module for our AMD machine retails for about $100 USD. 

Evolutions of BBH that start at larger separations would in principle require longer 
running times and more computer memory. Estimates of how much longer it would take 
to run extra orbits on any of our example runs can be obtained from the information 
given in the last column of Table [3j The memory requirements, on the other hand, 
will depend on the grid structure to be used and obviously the location of the outer 
boundaries. The memory requirements for runs such as ours would not change as long 
as the separation distance is not increased to more than 9.5M, which would yield an 
initial orbital period of more than 200M, thus leading to several more orbits before 
merger. The runs described in Table [2] were performed using quadrant symmetry. 
However, generic BBH with arbitrary masses and spins have to be simulated in full 
grids which, all things kept equal, would require about four times more memory and 
execution time. These requirements, however, can be reduced by adopting a different 
grid geometry. In Tichy and Marronetti [42] , generic BBH runs were performed using a 
grid Xr?=2[5 x 48 : 5 x 54] [hg = M/56.9 : OB = 240 M]. The performance details of one 
of these runs (completed on the Dual Intel Workstation) have been added in the last 
row of Table [3] where we see that, while slower than the previous runs, this simulation 
only takes a couple of days of execution time per orbit. 

Finally we performed the following comparisons. Firstly, we performed run C5 
also on a supercomputer (Cray XT3 MPP system at the Pittsburgh Supercomputer 
Center) using 32 processors. The execution on this machine was about 7 times faster 
(20 M/hr) than on our AMD Dual Opteron workstation. Secondly, we have also evolved 
model C5 of Table [2] with the LEAN code [37] . Because quadrant symmetry is not 
implemented in the current version of the grid driver of the LEAN code, this simulation 
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Table 3. Typical performance of the BAM code on a AMD Dual Opteron 2.2GHz 
workstation (runs PA, CI, C3 and C5,) and an Intel Dual Xeon 2.6GHz workstation 
(runs PI, P2, P3, P5 and C4). While these workstations are by themselves not 
particularly expensive, they were fitted with 8Gb (Opteron) and 16Gb (Xeon) of 
memory that increased their cost significantly. The last column shows the time (in 
days) it takes to evolve for one orbit using a nominal orbital length of 114M, obtained 
from the initial data set angular velocity. Due to the loss of the files, P6 performance 
could not be estimated. LC5 corresponds to a simulation identical to C5 but performed 
using the LEAN code [37] • The last row (FG) corresponds to a simulation using full 
grid presented in [42] , 

was performed using equatorial symmetry, but using four processors. The results are 
listed as LC5 in Table [21 We emphasize a few caveats in this last comparison. For 
instance, it was not possible to use identical grid setups due to the different types of 
symmetries used and the use of cell-centered and vertex-centered grids in BAM and 
LEAN respectively. Furthermore, the codes continue to undergo further development 
with likely improvements, in particular in the case of memory usage in BAM. Finally, 
the results are likely to be affected by the inclusion of further diagnostic tools, such as 
horizon finding. However, BAM and LEAN's performance are similar, both in terms 
of memory usage and speed, indicating that the critical aspect of these codes efficiency 
is the particular grid structure, more than in intrinsic coding details of the evolution 
equations. 

4. Code tests 

4-1. Convergence tests 

Table H] indicates the values obtained for the mass and angular momentum at the time 
the simulations were stopped (t F ). We also present the time of merger, estimated as 
the moment when the lapse function drops below a threshold value of 0.3. The error 
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Table 4. Values of the mass My and angular momentum Jy at the time the 
simulations from Tabled] were stopped (tp). The merger time {tu) is estimated from 
the time when the minimum value of the lapse drops below 0.3. The error bars are 
estimated from the variation in the last 100A/. Due to the loss of the files, P6 merger 
time could not be estimated. 

bars simply present the change of each quantity in the last 100M of the simulation (or a 
nominal value of 0.001, if the change is smaller than this threshold). The C-runs show 
larger changes in Jy than the P-runs due to the larger grid size used in the former 
simulations which requires longer evolution times for the volume integrals to settle. 

Figure [2] shows a comparison of My and Jy for different positions of the inner 
surface for run C5. The length of the inner cube side (d) is varied in multiples of 
d = 0.79M. Note for reference that the coordinate radius of the final black hole is 
~ 0.73M. The curves obtained for the smallest cube sizes show noise during the pre- 
merger stages which disappears when increasing the cube side. The spikes present at 
this stage are a numerical artifact caused by the crossing of the moving boxes of the 
x — z symmetry plane. After the merger, the curves settle to values that agree within a 
relative error of less that 0.2% for My and 0.04% for Jy. This seems to indicate that, 
when measuring the characteristics of the final black hole, the position of the inner cube 
does not affect greatly the results. 

Figure [3] compares the same quantities for two runs (PI and P5) that only differ 
in the size of the numerical grid. By adding two extra refinement levels to the outside 
of the grid used for PI, we moved the outer boundaries out by a factor of four. The 
main difference between these P-runs is the way the curves relax to their final values 
which agree to within a relative error of 0.01% (0.1%) for My (Jy). Again, this seems to 
show that the size of the grid does not affect significantly the final values of the global 
quantities. 

Figure H] shows a comparison of runs where the maximum grid resolution has been 
varied, leaving the rest of the grid characteristics identical. Note that the Jy curves 
present a downward trend that persists well after the merger. This trend gets less 
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Figure 2. Comparison of My and Jy for different inner surfaces. The curves 
correspond to run C5. The side of the inner surface cube is denoted by d. 
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Figure 3. Comparison of My and J v for two runs with different grid sizes. The 
curves correspond to runs PI (solid) and P5 (dashed). P5 is identical to PI, except 
for the presence of two additional low resolution outer levels which extend the outer 
boundaries by a factor of four. 
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Figure 4. Comparison of My and Jy for runs with different grid resolutions. The 
curves correspond to runs CI to C5. The values of My and Jy used for this test 
correspond to inner surface cubes of size 6.33M. 
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Figure 5. Convergence of My and Jy with grid resolution. The curves correspond 
to runs C2 to C5. The values of My and Jy used for these plots correspond to inner 
surface cubes of size 6.33M. 



pronounced for higher resolution and might be related to a similar coordinate drift 
observed in [12] . Fig. [5] presents convergence plots for the runs (72 to C5. The runs are 
grouped in two sets (72, (74, and (75 (thick upper curves) and (73, (74, and (75 (thin 
lower curves) and their differences are compared and scaled according to a putative 
fourth order convergence. The expected fourth order convergence is approached better 
by the set with higher resolution. 

We also compare two runs with different characteristics but identical maximum grid 
resolution. The curves from Fig. [6] correspond to runs PI (solid) and (73 (dashed). PI 



Binary black holes on a budget 



15 



0.98 
0.96 

> 

0.94 
0.92 
1.2 

> 1 
0.8 



"'"0 200 400 600 800 1000 1200 1400 

t/M 

Figure 6. Comparison of My and Jy from two runs using the dynamical variables 
4> (PI) and X (C3). 

was performed using <fi as the dynamical variable, Eq. ([1]) with 77 = 1.0/M in the recipe 
for the shift and it had outer boundaries placed at 244M. C3 used % as the dynamical 
variable, Eq. (j2J) with 77 = 2.0/M and outer boundaries at 727 M. The values of My 
and Jy used for these plots correspond to inner surface cubes of size 6.33M. As in Fig. 
El the difference in relaxation is mostly due to the difference in grid sizes. The relative 
difference between the global quantities at the end of these runs is 0.1% for My and 1% 
for Jy. 

Finally, in order to determine the values of the mass and angular momentum for 
the C-runs of Table EJ we used Richardson extrapolation with a polynomial of the type 

P(ho) = A + A x * h% + A 2 * h 5 . (12) 

Note that this formula assumes fourth order convergence which, according to Fig. El 
is only approximate for the runs of this paper. The extrapolation is performed at 
t = 1500M from runs C3, C4, and C5, giving My = 0.909 and Jy = 0.753. These 
values agree to within 3% of those reported in [T5] . 

4-2. Comparison of My and Jy with alternative estimates of the mass and angular 
momentum of the final black hole 

Alternative estimations of the mass and angular momentum of the final black hole 
can be derived from evaluating the surface integrals (0 ED at finite radii. Figure [7] 
compares the values of My and Jy with integrations performed on spherical surfaces at 
coordinate radii 30M, 50M and 100M. The mass curves corresponding to the spherical 
surface integrations have very large errors that behave in a non-systematic way when 
the extraction radius is increased. At the same time, the calculations of the angular 
momentum using spherical surface integrals, while noisy before the merger, agree with 
Jy to within a 0.001% relative error at the end of the simulation. 
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Figure 7. Comparison of My and Jy (solid) and spherical surface integrations at 
radii 30M (dotted), 50M (dashed), and 100M (dashed-dotted). The curves correspond 
to run C4. 

Figure [8] shows the radiated energy calculated from the flux through spherical 
surfaces using the Newman-Penrose scalar ^4 (see Eq. (52) in [12]) at radii 15M, 
30M, 60M, 80M for run C4. From the radiated energy, we generated an estimate for 
the time evolution of the mass by subtracting those curves from the initial M ADM as 
reported in Table [1] (dashed curves). These curves show a downward drift that is more 
pronounced with the proximity of the spherical surface to the center of the grid. The 
result obtained for radius 80M is the one that appears to be the least affected by this 
effect and it agrees with My within a relative error of about 1%. 

The Christodoulou formula, valid for stationary axisymmetric spacetimes like the 
Kerr geometry of the final black hole, gives the following relation between the angular 
momentum J, mass M and irreducible mass of a black hole M irr 



J = 2M irr ^JM 2 - Mf„ . (13) 

Here 



Mi, 



.4 



H 



V 167T 

is determined from the proper area Ah of the apparent horizon (for the run C4 
Mi rr ~ 0.825M). We verified relation ( TL31) for several of the runs of Table [2] and found 
it to be satisfied to a relative error of less than 0.2%. 

5. Results and conclusions 



The main goal of this paper is to show the ability of the BAM code to perform certain 
simulations of binary black holes with relatively modest computer resources. The 
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Figure 8. Comparison of My (black solid) and mass estimates from the energy 
radiated through spherical surfaces at different radii (color dashed) . The latter curves 
were generated by subtracting the radiated energy from the initial M ADM . The curves 
correspond to run C4 and inner surface cube size 6.33Af. 



simulations presented here were performed on dual processor workstations that have 
been outfitted with at least 8Gb of memory. Machines like these retail at the moment 
of publication for less than $3000 USD, making them easily accessible to any research 
group. 

Our runs were based on the same initial state: one corresponding to two identical 
black holes with intrinsic spin parameters S/m 2 = 0.75 and spins parallel to the orbital 
angular momentum and which started out at a coordinate separation of 6M. This data 
set is similar to the one used by Campanelli et al. fTS]. Our results for the mass and 
angular momentum of the final black hole agree to within 3% of those of [15J. Better 
accuracy could be achieved with higher resolution runs, however they would also demand 
larger computer resources. We are currently testing different grid structures (i.e., varying 
moving boxes sizes, number, etc.) that improve these runs accuracy without increasing 
the computational burden (see, for instance [12] )• We performed high spatial resolution 
simulations (of up to M/160) and large grid size (up to outer boundaries at 975M) to 
test the moving punctures method's stability and robustness. None of the runs showed 
any signs of exponentially growing instabilities; they were stopped due to the long real- 
time duration of the simulations. In one of our test cases, the simulation was run for 
the equivalent of more than 20 orbital periods. 

The simulations discussed in this paper were performed using quadrant symmetry 
which require about four times less memory and computer time than non-symmetric 
scenarios. An improved choice of grid structure can (to some extent) minimize these 
requirements (see [12], for non-symmetric simulations using workstations). BAM 
capabilities for efficient use of computer resources permit the exploration of BBH 
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parameter space by enabling low resolution simulations that only require workstations 
or one or two nodes per run in local Bewoulf clusters, where many of such runs can 
be done simultaneously. Ongoing code optimization is currently enhancing the code 
performance with regard to memory and cpu usage. Nevertheless, it is clear that today's 
most demanding binary simulations still require supercomputer resources. However, 
given the rapid growth of computer power and high efficiency codes such as BAM, even 
these simulations might be within the reach of workstation resources in the next couple 
of years. 
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